clear;
# constantes globales
N = 1000;
clases = 7;
seed = 10;
#seed = 1105;
#seed = 31416;

for i=1:clases
	n(i) = i;
	O(i) = 0;
	O2(i) = 0;
	ON(i) = 0;
	X(i) = 0;
	F(i) = 0;
	Fn(i) = 0;
	Fn2(i) = 0;
endfor

x(1) = seed;
for i=1:2 * N + 1
	x(i+1) = mod(25214903917 * x(i) + 11, 2^48);
endfor

for i=1:N
	u(i) = (2^27 * floor(x(2*i) / 2^22) + (floor(x(2*i+1) / 2^21))) / 2^53;
	Norm(i) = (u(i)^0.135 - (1 - u(i))^0.135) / 0.1975;
	for j=1:clases
		if (u(i) > (1 / clases) * (j-1) && u(i) < (1 / clases) * (j))
			O(j) = O(j) + 1;
		endif
	endfor
	for j=1:clases
		if (Norm(i) > (1 / clases) * (j-1) && Norm(i) < (1 / clases) * (j))
			ON(j) = ON(j) + 1;
		endif
	endfor
endfor

for i=1:N-1
	X1(i) = u(i);
	Y1(i) = u(i+1);
endfor

for i=1:N-2
	X2(i) = u(i);
	Y2(i) = u(i+1);
	Z2(i) = u(i+2);
endfor

# TEST CHI CUADRADO
SumChi = 0;
printf("TEST CHI CUADRADO\n");
printf("Clases\tO(i)\tE(i)\t\tO(i)-E(i)\t(O(i)-E(i))^2\t\t(O(i)-E(i))^2 / E(i)\n");
for j=1:clases
	E(j) = N / clases;
	DifCuad(j) = ((O(j)-E(j))^2);
	Chi(j) = DifCuad(j) / E(j);
	SumChi = SumChi + Chi(j);
	printf("%d\t%g\t%g\t\t%g\t\t%g\t\t\t%g\n",j,O(j),E(j),O(j)-E(j),DifCuad(j),Chi(j));
endfor
printf("\nEl estadistico Chi Cuadrado es %g\n\n",SumChi);

# TEST KS
u2 = sort(u);
for i=1:N
	for j=1:clases
		if (u2(i) > (1 / clases) * (j-1) && u2(i) < (1 / clases) * (j))
			O2(j) = O2(j) + 1;
		endif
	endfor
endfor
Dpos = 0;
Dneg = 0;
D = 0;
printf("TEST KS\n")
printf("Xi\ti/n\ti/n-Xi\tXi-(i-1)\n");
for j=1:clases
	X(j) = O2(j) / N;
	F(j) = (1/clases) * j;
	Fn(j) = F(j) - X(j);
	if (j == 1)
		Fn2(j) = X(j);
	else
		Fn2(j) = X(j) - X(j-1);
	endif
	if (Fn(j) > Dpos)
		Dpos = Fn(j);
	endif
	if (Fn2(j) > Dneg)
		Dneg = Fn2(j);
	endif
	printf("%g\t%g\t%g\t%g\n",X(j),F(j),Fn(j),Fn2(j));
endfor
D = Dpos;
if (Dneg > D)
	D = Dneg;
endif
printf("\nEl estadistico Dpos es %g y el Dneg es %g, entonces el KS es %g\n",Dpos,Dneg,D);

# SIMULACION DE MONTECARLO
ENorm = 0;
for i=1:N
	ENorm = ENorm + abs(Norm(i));
endfor
ENorm = ENorm / N;
printf("\nSIMULACION DE MONTECARLO\n");
printf("Si Z~N[0,1], la E(Z) = %g\n", ENorm);

# PLOT HISTOGRAMA UNIFORME
#hist(u,clases);
#ylabel("Cantidad");
#xlabel("Clases");
#print -deps "~/Documents/ITBA/ss/tp4/HistUniforme.eps";

# PLOT HISTOGRAMA NORMAL
#hist(Norm,clases+4);
#ylabel("Cantidad");
#xlabel("Clases");
#print -deps "~/Documents/ITBA/ss/tp4/HistNormal.eps";

# PLOT PARES
plot(X1,Y1,".");
#xlabel("Ui");
#ylabel("Ui+1");
#print -deps "~/Documents/ITBA/ss/tp4/Pares.eps";

# PLOT TERNAS
#plot3(X2,Y2,Z2,".");
#xlabel("Ui");
#ylabel("Ui+1");
#zlabel("Ui+2");
#print -deps "~/Documents/ITBA/ss/tp4/Ternas.eps";
